#ifndef EST_POLYA_FIXEDPOINT_CC
#define EST_POLYA_FIXEDPOINT_CC

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include <gsl/gsl_sf_psi.h>
#include "def.hpp"
#include "utils.hpp"
#include "io.hpp"

void est_polya_fixedpoint(int**& data, int*& n, double*& estimated_alpha, int& K, int& N){
  double* a_buffer = new double[K];
  double b_buffer = 0.0;
  double* alpha_buffer = new double[K];
  double* temp_buffer = new double[K];

  for(int k=0; k < K; k++){
    alpha_buffer[k] = estimated_alpha[k];
  }

  for(int iter =0; iter < 1000; iter ++){
    // some preparation
    double temp_sum_alpha = get_sum(alpha_buffer,K);
    for(int k=0; k < K; k++){
      a_buffer[k] = 0.0;
      temp_buffer[k] = 0.0;
    }
    b_buffer = 0.0;

    for(int i =0; i < N ; i++){
      for(int k=0; k < K; k++){
	a_buffer[k] = a_buffer[k] + gsl_sf_psi(data[i][k] + alpha_buffer[k]);
      }
      b_buffer = b_buffer + gsl_sf_psi(n[i] + temp_sum_alpha);
    }

    for(int k=0; k < K ; k++){
      temp_buffer[k] = alpha_buffer[k] * ((a_buffer[k] - N * gsl_sf_psi(alpha_buffer[k])) / (b_buffer - N * gsl_sf_psi(temp_sum_alpha)));
    }

    // l1 distance
    double error = l1_distance(temp_buffer, alpha_buffer,K);
    if(error < 1e-8)
      break;
    for(int k=0; k < K ; k++){
      alpha_buffer[k] = temp_buffer[k];
    }

  }

  for(int k=0; k < K ; k++){
    estimated_alpha[k] = alpha_buffer[k];
  }

  delete[] a_buffer;
  delete[] alpha_buffer;
}

#endif
